Linear scaling and noisy forces

This paper by Krajewski and Parrinello,
Linear scaling electronic structure calculations and accurate sampling with noisy forces, has two intriguing ideas.

They compute electronic structure using a new linear scaling method that involves sampling, and thus results in noisy forces. Those forces are used to drive a molecular dynamics simulation of the nuclei (note that one of the authors is Parrinello).

For the electronic structure part of the problem, they replace a determinant in the independent particle grand potential with an integral over some fields. These fields can be sampled using a Langevin equation. This Langevin equation apparently involves only matrix-vector multiplication with a sparse matrix, which delivers the scaling.

The second part (starting in the paragraph before equation 11) is how to deal with the noisy forces in molecular dynamics. Their solution for this is to again use a Langevin equation, and to choose the size of the damping terms to cancel out the noise on average. In this way they obtain correct sampling of the Boltzmann distribution.

Polarization

Calculating observable properties of systems has been a weak point of QMC methods. The paper, Dielectric Response of Periodic Systems from Quantum Monte Carlo Calculations by Umari, Williamson, Galli, and Marzari remedies some of that deficiency by demonstrating a calculation of polarization using DMC.

"Programming in Math" updated

The paper on Programming in Mathematical Notation has been updated. (see previous blog entry about Programming in Math)

I've added a page about integration by random sampling. This example adds Greek characters, and a foreign function interface for the random number generator. Also, the translation of several separate sums (defining average, variance, etc) into C++ requires the sums to be merged into a single loop.

[Edit 2/24/2007 - update link]

Cluster algorithm

In the paper, Cluster algorithm for pairwise-interacting particles, Stephen Whitelam and Phillip Geissler describe an algorithm for creating a cluster of particles to move as a group.

The cluster starts with a randomly chosen particle and grows by adding particles with some probability (based on the Boltzmann factor). Particles that strongly interact with cluster are more likely to be part of the cluster. Then a normal trial move is made of the cluster as a whole.

Importantly, the potential and temperature chosen for this operation are not the same as the potential and temperature of the system. In fact, the cluster size can be controlled by adjusting the temperature (at infinite temperature each cluster is only a single particle). This cluster temperature is varied during the simulation run to sample moves with different cluster sizes.

See the paper for more details about controlling the range of energies involved in the cluster (and the cluster interface energies).

Programming in Math

In a previous post on Solving science problems, I commented on the process of scientific computation. One of the steps was converting the solvable equations into a computer program.

As part of making that step more transparent, I've been working on method for Programming in Mathematical Notation (works best with MathML capable browser).

The program takes a content MathML file as input and outputs a presentation MathML file for nicely rendered display in the browser, and a C++ file to compile and run.

[Edit 2/24/2007 - update link to example]

Fast variance optimization

VMC variance optimization usually uses the reweighted variance.

What is reweighting? Consider a set of particle positions obtained from sampling a wavefunction at a given set of parameters. The energy (and the variance of the energy) explicitly depends on the wavefunction parameters by the local energy formula, and implicitly through the particle positions determined by sampling the wavefunction. When optimizing, one would like to vary the parameters without going through all the trouble and expense of sampling new positions. Computing the energy (or variance) change due to the local energy piece is easy - just put the new parameters in and recalculate. Reweighting is the process of correcting the implicit dependence of the energy (or variance) on the wavefunction parameters.


N. D. Drummond and R. J. Needs describe a VMC optimization method in their paper, A variance-minimization scheme for optimizing Jastrow factors, that uses only the explicit dependence of the variance on the parameters (ie, they use the un-reweighted variance). They restrict their scheme to linear parameters in the Jastrow factor, and this makes the resulting parameter surface a quartic. Also, this scheme only requires a single sum over electron positions (the requisite factors can be accumulated during a single run), making the optimization step very fast.

Even though this scheme uses a less accurate approximation for the variance, it appears to give very good results. In fact, the resulting energies for a neon atom are better using the unreweighted variance for optimization than from using the more traditional reweighted variance.

Solving science problems

What is the ideal programming language for MC simulations?

What if we could design a language expressly for MC and QMC simulations - what would it look like? As a first step, let's put programming in the context of solving the whole problem.

Steps to solve a computational science problem

1. Define the problem (ask the question)
   
2. Write down the equations modeling the system under investigation
   
3. Using mathematical transformations and approximations, reduce the equations to a solvable form (numerical analysis)
   
4. Convert into a computer program and run it
   
5. Undertake verification (are steps 2,3, and 4 correct) and validation (are step 1 and the approximations in step 2 correct)
   

Developing techniques for managing each of the steps would be a great help. Current programming languages only deal with step 4. Computer algebra systems (Mathematica, Maple, Maxima, etc) can help with steps 2 and 3, but don't really help in making the jump to a large scale program.

Currently the transformations from steps 2 through 4 are done by hand. I know I make many mistakes in this process. It would be a boon to my productivity if the computer could verify that all the minus signs and factors of 2 are correct.

The verification step often involves understanding the effects of various approximations in step 3. If the program is described as a series of transformations on the equations, hopefully changing the approximation would be a simple change to the equations, and the rest of the process would automatically generate a correct program. This would make it easier to test approximations.

This doesn't say much about the ideal programming language, but it does give some idea of the problem to be solved. I'll have more to say about the ideal programming language later.

Programming Language Notes

On the Sun website, there's an interview with Guy Steele. He talks a little bit about Fortress, which tries to make programs looks more the original problem domain by allowing more mathematical notation. I will be interested to see the results.

Guy Steele is also the principal investigator for the Programming Language Research group. On that page there is a paper titled Object-Oriented Units of Measurement. It's a look the issues involved in making units work in a Java-like OO language.

The last note is an ACM Queue article by Gregory Wilson about Extensible Programming. The article points out that there are advantages to storing the program text in a more structured form, including being able to include diagrams or math notation in the source code.

Whack A Mole

In QMC, a central role is given to the local energy, defined as (1/psi)H psi, where psi is the trial wave function. If psi is an eigenstate, then H psi = E psi by definition, and the local energy becomes a constant. Just as lower energy trial functions are deemed better than higher energy ones, those with fewer fluctutations are better than those with more fluctutations.

One could imagine looking at the fluctuations in the local energy, and using those to determine where best to adjust the trial wave function to reduce the fluctutations (side note: one would want to use psi^2 times the local energy, since fluctuations are weighted by psi^2 in the integrals). This seems similar to the game of Whack-A-Mole, where one attempts to hit mechanical moles that pop up randomly for a brief period of time. Maybe the method could be called Whack-A-Local-Energy (WHALE).

Unfortunately for my naming scheme (and fortunately for everyone else), the method described already exists, and it's called the Energy Fluctutation Potential (EFP) method. It was first described by Stephen Fahy in "Quantum Monte Carlo Methods in Physics and Chemistry" [NATO ASI Ser. C 525 101, 1999; M. P. Nightingale and C. J. Umrigar, eds]. Further papers describing the method:

• Optimal orbitals from energy fluctuations in correlated wave functions, by Fahy and Filippi
  
• Optimization of inhomogeneous electron correlation factors in periodic solids by Prendergast, Bevan, and Fahy
  
• Optimized Jastrow-Slater wave functions for ground and excited states: Application to the lowest states of ethene by Schautz and Filippi
  

The EFP method involves an iteration that starts with a one-body solver (HF or LDA). The resulting orbitals are used in a VMC calculation, and the fluctuations in the local energy form a extra potential that is fed back into the one-body solver to get new orbitals, and the iteration repeats until a self-consistent solution is reached.

(Name note: this paper by Umrigar and Filippi calls it the Effective Fluctuation Potential method. This paper also gives another reference I didn't put on the above list)

Store, manage, and share your reading list

The cleverly-named del.icio.us site describes itself as "social bookmarks." The "bookmarks" part allows users to keep track of links online, include a short description, and add an arbitrary number of one-word tags to help categorize the link. With the
"social" part, people's bookmarks are available to everone else.


Here's my list of bookmarks. (I just started, it's not very large)

For the academic world, Richard Cameron created CiteULike. It allows assigning category tags to papers, and has the social aspect of sharing links to papers. Perhaps most importantly, it can automatically import bibliographic information from supported publishers (including arXiv) and export the information to other formats (BibTex, EndNote).


Here's my list of articles. (Once again, very small so far)

The link to CiteULike was found via Corante

I don't know whether either of these will be a useful tool for research, but it seems worthwhile to try them out. When I started grad school in 1993, I went to the library to photocopy articles. By the time I finished in 2000, I went to the library website to download articles. That seems like a large shift, but I suspect it's not just because it makes it faster and more efficient to get the articles - it's the linking, management, and annotations represented by these sorts of applications that will make a bigger change in the future.

MathMath

A wiki seems like a good way to organize information. Blogs work nice for recording thoughts and links in chronological order, but aren't so good for organizing ideas topically or for revising items. Proto-papers or latex files with background information have been used. The easily edited environment of a wiki seems like it would function well, even if were only used personally (part of the wiki idea is open collaborative editing). The wiki text markup is much lighter than HTML, so it would be easier to deal with than maintaining a set of web pages.

The real issue is integrating mathematics. Mozilla and derivatives can handle MathML (and IE can with a plugin), so that seems like a good route (inserting pictures is marginally okay for converting papers to HTML, but is not as good for a more interactive environment). The only problem with MathML is the rapid onset of RSI if you actually write equations with it.

There are a couple of solutions, both convert TeX or vaguely TeX-like input to MathML. They are itex2mml and AsciiMathML. The second solution seems rather nice, as it uses javascript to convert the TeX-like notation to MathML when the page loads into the browser.

For a wiki, I've been looking into MoinMoin. At first I tried hacking itex2mml into it, and it took quite a few changes to make everything proper XHTML. Then I found the page that describes integrating AsciiMathML into MoinMoin, and that was much smoother.

AsciiMathML also has a nice way of adding new symbols, that doesn't require changing the AsciiMathML.js file.
In keeping with title of this weblog, the Unicode number for hbar is 0x210F. (although it usually only appears once in a paper and is promptly set to "1") The code for adding it is

<script type="text/javascript">
AMsymbols = AMsymbols.concat([
{input:"hbar", tag:"mo", output:"\u210F", tex:"hbar"}
]);
</script>

Chemistry Textbooks

Jack Simons has a couple of introductory chemistry textbooks on his website. They are An Introduction to Theoretical Chemistry and Quantum Mechanics in Chemistry (with Jeff Nichols).

Stochastic Matrix Inversion

I ran across this paper, A Monte Carlo algorithm for efficient large matrices inversion. There are a couple of techniques presented, one which is similar to the iterative Gauss-Seidel algorithm.

QMC methods often use a matrix to enforce antisymmetry conditions on a fermionic wavefunction. The determinant is needed for the wavefunction itself, and the inverse matrix is used for the local energy.

I wonder if a stochastic matrix algorithm would be useful in the QMC calculations. My guess is that it would not be (and if there is an efficiency gain, it might be offset by the added complexity), but it would be interesting to try, if for no other reason than ideological purity (if you're going to use a stochastic method, might as well make it as stochastic as possible!)

Other interesting questions arise. The locations of the all-important wavefunction nodes will now be fuzzy. What consequences does this have for a fixed-node algorithm?

Inverse matrices also appear in Sorella's Stochastic Reconfiguration method for optimization. Another use?

On a side note, I tried looking for background on Google, but searching for "stochastic matrix inversion" returned results relevant to inverting a stochastic matrix rather than stochastic methods for inverting a matrix.

Nifty math and physics applets

Paul Falstad has a page of math and physics applets. (found via Alan Aspuru-Guzik's weblog)